# local paths ----
# inputs
ships_shp <- here("data/raw/ships/boem_gom_grid_pl_albers/boem_gom_grid_pl_albers.shp")
ships_xls <- here("data/raw/ships/BOEM GoMex 2020 BiOp AIS Data for Rices Whale Risk Analysis.xlsx")
whales_shp <- here("data/raw/whales/Rices_Whale_Monthly_Density.shp")
# outputs
study_geo <- here("data/study.geojson")
units_geo <- here("data/units.geojson")
ships_geo <- here("data/ships.geojson")
whales_geo <- here("data/whales.geojson")
ships_csv <- here("data/ships.csv")
ships_avg_csv <- here("data/ships_avg.csv")
whales_csv <- here("data/whales.csv")
# read data ----
# transform all to geographic projection (4326)
if (!all(file.exists(c(
study_geo, units_geo, ships_csv, whales_csv)))){
# ships
grd_ships <- read_sf(ships_shp) |>
st_transform(4326)
# whales
hex_whales <- read_sf(whales_shp) |>
st_transform(4326)
# study area ----
# as intersection of all dissolved (st_union) inputs: whales, ships, eez
if (!file.exists(study_geo)){
# eez
us_eez <- mr_features_get(
type = "MarineRegions:eez",
featureID = "eez.281") |>
geojson_sf() |>
st_transform(4326) # mapview(us_eez)
# intersect all
ply_study <- st_union(grd_ships) |>
st_intersection(
st_union(hex_whales)) |>
st_intersection(
us_eez)
# write geojson
write_sf(ply_study, study_geo)
}
ply_study <- read_sf(study_geo)
# units: intersection of ships, whales, eez ----
if (!file.exists(units_geo)){
grd_ships |> write_sf(here("data/ships.gpkg"))
hex_whales |> write_sf(here("data/whales.gpkg"))
# ran intersection in QGIS
stop("run intersection in QGIS manually")
# processing.run("native:intersection", {'INPUT':QgsProcessingFeatureSourceDefinition('/Users/bbest/Github/ecoquants/ricei/data/whales.gpkg|layername=whales', selectedFeaturesOnly=False, featureLimit=-1, flags=QgsProcessingFeatureSourceDefinition.FlagOverrideDefaultGeometryCheck, geometryCheck=QgsFeatureRequest.GeometryNoCheck),'OVERLAY':QgsProcessingFeatureSourceDefinition('/Users/bbest/Github/ecoquants/ricei/data/ships.gpkg|layername=ships', selectedFeaturesOnly=False, featureLimit=-1, flags=QgsProcessingFeatureSourceDefinition.FlagOverrideDefaultGeometryCheck, geometryCheck=QgsFeatureRequest.GeometryNoCheck),'INPUT_FIELDS':[],'OVERLAY_FIELDS':[],'OVERLAY_FIELDS_PREFIX':'','OUTPUT':'ogr:dbname=\'/Users/bbest/Github/ecoquants/ricei/data/ships_whales.gpkg\' table="ships_whales" (geom)','GRID_SIZE':None})
ply_units <- here("data/ships_whales.gpkg") |>
read_sf("ships_whales") |>
st_intersection(ply_study) |>
select(
hexid = HEXID,
grdid = grid_id) |>
st_make_valid() |>
mutate(
area_m2 = st_area(geom))
write_sf(ply_units, units_geo)
file_delete(dir_ls(here("data"), glob = "*.gpkg"))
ply_ships <- ply_units |>
group_by(grdid) |>
summarize()
write_sf(ply_ships, ships_geo)
ply_whales <- ply_units |>
group_by(hexid) |>
summarize()
write_sf(ply_whales, whales_geo)
}
ply_units <- read_sf(units_geo)
# tables (*.csv): ships, whales limited to study ----
if (!file.exists(ships_csv)){
tbl_ships <- read_excel(ships_xls) |> # 313,200
rename(grdid = grid_id) |>
filter(
grdid %in% unique(ply_units$grdid))
write_csv(tbl_ships, ships_csv)
tbl_ships_avg <- tbl_ships |>
group_by(grdid) |>
summarize_all(sum) |>
select(-yr, -mo)
write_csv(tbl_ships_avg, ships_avg_csv)
}
if (!file.exists(whales_csv)){
tbl_whales <- hex_whales |>
st_drop_geometry() |>
clean_names() |>
filter(
hexid %in% unique(ply_units$hexid)) |>
mutate(across(where(is.numeric), ~na_if(., -9999))) |>
filter(if_any(!hexid, ~ !is.na(.))) |>
rowwise() |>
mutate(
avg_n = mean(c_across(ends_with("_n")), na.rm = T))
write_csv(tbl_whales, whales_csv)
}
}
ply_study <- read_sf(study_geo)
ply_units <- read_sf(units_geo)
ply_ships <- read_sf(ships_geo)
ply_whales <- read_sf(whales_geo)
tbl_ships <- read_csv(ships_csv)
tbl_ships_avg <- read_csv(ships_avg_csv)
tbl_whales <- read_csv(whales_csv)